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FOREWORD 


The  following  articlee  wer^j  translated  becauae  they  are  relevant 
to  computatlona  performed  for  the  Navy  Space  Surveillance  Project  and 
the  Transit  Navigation  Satellite  Project.  The  articlee  should  also 
be  of  interest  to  others  concerned  with  the  subject  of  errors  in  numeri¬ 
cal  integration. 
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THE  ESTIMATION  OF  THE  ERROR  RESULTING  FROM  NUMERICAL  INTEGRATION  OF 
THE  EQUATIONS  OF  CELESTIAL  MECHANICS 

V.  F.  Myachin 

This  article  presents  the  result  of  the  application  of  a  general 
theory  of  estimation  of  the  numerical  integration  error  of  differen¬ 
tial  equations,  proposed  by  Professor  S.  M.  Lozinsky,  to  equations 
of  undisturbed  motion  of  celestial  mechanics.  In  applying  the  theory, 
consideration  was  given  to  the  random  character  of  the  round-off 
error.  Here,  in  particular,  it  is  shown  that  the  new  theory  gives 
qualitative  confirmation  of  the  well-known  result  of  Brouwer,  which 
asserts  that  after  a  sufficiently  large  number,  k  ,  of  integration 
steps  the  error  in  the  coordinates  of  elliptic  motion  has  the  order 

_3_ 

of  growth  Jt* . 


Introduction 

Let  us  consider  a  system  of  equations  of  disturbed  motion: 

i  =  /?„ 

(A) 

As  is  generally  known,  the  most  prevalent  method  of  approximating 
the  solution  of  such  equations  is  by  numerical  integration;  in  addi¬ 
tion  it  is  known  that  the  application  of  any  formula  of  numerical 
integration  inevitably  leads  to  an  accumulation  of  error  in  the  com¬ 
puted  approximate  solution,  which  arises  in  two  ways: 

1)  use  in  the  integration  formulas  of  a  finite  number  of  terms; 

2)  errors  in  initial  conditions  and  round-off  in  each  step  of 
integration. 
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At  the  present  time  we  do  not  have  available  strict,  practically 
acceptable  estimates  of  the  error  of  numerical  integration  of  differ¬ 
ential  equations;  however,  there  exist  so-called  rough  estimates, 
which  were  obtained  by  rough  reasoning. 

Thus,  Newcomb  (1898)  regarded  the  error  in  the  numerical  integra¬ 
tion  of  any  differential  equation  as  the  result  of  the  sum  of  the 
round-off  errors  at  each  step  (by  analogy  with  the  error  of  appi'oxi- 
mate  computation  of  a  definite  integral);  counting  these  errors  as 
independent  random  variables  (in  the  sense  of  the  theory  of  proba¬ 
bility),  he  came  to  the  conclusion  that  the  accumulation  of^error 

after  k  steps  of  numerical  Integration  is  of  the  order  of  k*  for 

2. 

equations  of  the  second  order  and  k*  for  equations  of  the  first  order. 

Applying  this  assertion  to  the  estimation  of  the  accumulation  of 
round-off  errors  in  the  method  of  special  perturbations,  one  can  draw 
a  conclusion  about  the  advantage  of  obtaining  the  perturbations  in* 
the  elements  compared  with  their  computation  by  the  method  of  Cowell 
or  Encke  (i.e.,  in  the  coordinates).  In  the  first  case  the  perturba¬ 
tions  are  obtained  as  a  result  of  single  integration  for  all  elements 
except  the  mean  longitude,  which  requires  double  integration,  and, 
consequently,  only  in  this  element  is  there  a  serious  accumulation  of 
error.  If,  however,  the  perturbations  are  computed  by  the  method  of 
Cowell  or  Encke,  then  for  every  coordinate  there  will  be  an  equation 
of  the  second  order  and  the  expected  error  in  each  coordinate  will  be 
of  the  order  of  *  . 

Brouwer  (1937),  however,  showed  that  this  conclusion  is  incorrect 
and  that  for  elliptic  motion  the  accumulation  of  errors  in  the 

coordinates  is  equivalent  to  an  error  in  mean  longitude  proportional 

—  L 

to  k* ,  and  to  errors  in  the  other  elements  proportional  to  k*  . 
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Consequently,  for  circular  motion,  where  the  mean  motion  is  a  con¬ 
stant  quantity,  the  error  in  the  coordinates  by  the  integration 
method  of  Cowell  or  Encke  is  proportional  to  k*  .  In  this  connection 
Brouwer  first  noticed  the  essential  difference  between  the  approxi¬ 
mate  computation  of  a  definite  integral  and  the  numerical  integration 
of  differential  equations;  in  the  latter  case  the  computation  of  each 
tabular  value  uses  the  result  of  a  prior  integration,  so  that  the 
error  in  this  case  cannot  possibly  be  considered  as  the  result  of 
simple  summation  of  errors  at  each  step. 

Still,  because  of  loose  reasoning,  the  result  of  Brouwer  on  the 
quantitative  side  appears  rough:  in  particular,  it  does  not  reflect 
the  experimentally  known  fact  that  the  error  has  an  oscillatory 
character;  moreover,  at  a  sufficiently  large  number  of  steps  of 
integration  the  formulas  give  an  underestimate,  i.e.,  the  true  error 
exceeds  prediction. 

At  a  later  time  Professor  S.  M.  Lozinsky  proposed  a  number  of 
general  theorems  for  estimating  the  error  in  numerical  integration 
of  systems  of  differential  equations. 

The  present  paper  presents  the  result  of  the  application  to 
System  (A)  of  Lozinsky's  so-called  "linearized  estimate"^,  which  was 
somewhat  modified  and  adapted  in  the  interest  of  our  problem.  It  is 
shown  that  il',  following  Newcomb  and  Brouwer,  one  considers  the 
round-off  errors  as  independent  random  variables,  then  for  the  accumu¬ 
lation  of  error  in  the  coordinates  of  elliptic  motion  the  indicated 
estimate  gives  the  order  of  A*,  and  for  circular  motion  the  order 
of  , 

^At  the  present  time  this  estimate  is  in  press.  By  permission  of  the 
author  we  borrowed  it  from  a  course  of  lectures  "Approximate  Solu¬ 
tion  of  Differential  Equations",  read  by  S.  M.  Lozinsky  at  Leningrad 
University  in  1955-1957. 
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This  method  permits  one  to  take  account  not  only  of  the  influence 
of  round-off  errors  on  the  error  in  the  coordinates,  but  also  the  effect 
produced  by  errors  in  the  initial  conditions  and  by  the  discard  of  the 
remainder  term  in  the  integration  formvilas.  Moreover,  it  is  possible 
to  keep  track  of  the  variation  of  the  error  separately  in  each  coordi¬ 
nate. 

The  exposition  is  based  upon  an  example  of  Cowell's  method. 

1.  A  FEW  GENERAL  STATEMENTS  ABOUT  MATRICES  AND  SYSTEMS  OF  LINEAR 
DIFFERENTIAL  EQUATIONS 


Let  there  be  given  the  rectangular  matrix; 


.4  = 


Out  On, 


•»  fll» 
•*  fl2» 


Omit  Qm2»  •  •  •(  Omn 


We  shall  use  the  following  notation: 

A  matrix  consisting  of  one  column  is  called  a  vector.  Let  *  be 
a  vector  with  components  ...,  x^”}. 

We  shall  denote: 

=  . x<‘)=U)(0. 

By  |jt|  we  shall  denote  a  vector,  defined  by  the  equation 

Finally,  we  introduce  the  norm  of  the  vector  je,  denoted  by  |jt} 
and  defined  by  the  equation 

lje(  =  m"ax|Un. 


A  vector  will  be  considered  as  a  special  case  of  a  matrix. 

For  a  given  matrix  A  we  denote  by  A  the  matrix  defined  by  the 
equation 

{A),j^(A)j,i 


it  is  called  the  transpose. 

The  transpose  of  a  vector  is  a  matrix  consisting  of  one  row: 

.... 

A  matrix  having  the  same  number  of  rows  and  columns  is  called  a 
square  matrix.  The  elements  {y4)«  of  a  square  matrix  A  form  the 
principal  diagonal. 

Operations  upon  matrices: 

1.  Addition  of  two  matrices.  By  definition: 

2.  Multiplication  of  a  matrix  /I  by  a  number  X  .  By  definition: 

(X/l}iyS!X{^},y. 


3.  Matrix  multiplication.  Let  a  matrix  A  have  n  columns,  and  a 
matrix  B  have  n  rows;  then  by  definition: 

\AB),j=  l(Ah^{B)„j. 

If  the  number  of  columns  of  A  is  not  equal  to  the  number  of  rows 
of  B,  then  the  product  <45  is  not  defined.  Matrix  multiplication 
possesses  the  associative  property,  Siid  multiplication  and  addition 
possess  the  distributive  property: 

{AB)C=A{BC), 

C{A-^^B)D  =  CAD  CBD, 
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if  all  indicated  products  are  defined. 

The  operation  of  matrix  multiplication,  generally  speaking,  does 
not  have  the  commutative  property: 

AB^BA. 

For  example,  for  the  vector  ....  x*"))  and  its 

transpose  !t  we  have: 

I  .1.1 


but 


XX  = 


x0)\  x^*)xi‘), 


X<‘>X<") 


(1) 


4.  The  transpose  of  the  product  of  two  matrices: 

(A~B)  =  BA.  (2) 

5.  Matrix  differentiation.  Let  A{t)  be  a  matrix  representing  a 
function  of  the  scalar  argument  /;  then  by  definition: 


6. 


Integration. 


By  definition; 


\A{f)<U' 


\{A{f)),,df. 


U.  )o  '» 

Henceforth,  we  will  speak  only  of  square  matrices.  A  matrix  whose 
elements  are  all  zero  is  called  null  and  is  denoted  by  0.  A  matrix 
whose  elements  on  the  principal  diagonal  are  equal  to  unity,  and  all 
the  rest  are  zero,  is  called  a  unit  matrix  and  denoted  by  /. 


1. 

0.  . 

0 

/  = 

0, 

1, . 

..  0 

0. 

0.  . 

..  1 
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For*  any  matrix  A  it  in  true  that 

AI^IA  —  A. 

The  matrix  S,  possessing  the  property 

SS^I, 

Is  called  orthogonal. 

Consider  now  the  system  of  linear  differential  equations 

y(0—  2 

Introducing  the  vector 

.V  {«/".  !/■',  . . j/"*) 

and  the  matrix 

P(/) 

we  write  this  system  in  the  form 

rr-pm-  (3) 

We  introduce  the  matrices  of  order  n  V {i^  t),  defined 

by  the  conditions: 

0=:P{t)U.  U{t„t,)  =  I,  t/(/„O  =  0 

i^=p{t)v.  v{t„  g=o.  /,)=/, 

which  are  called  respectively  the  first  and  second  fundamental  solu¬ 
tions  of  System  (3) .  The  columns  of  these  matrices  represent  the 
linear  independent  solutions  of  System  (3). 
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Knowing  the  fundamental  solutions  of  System  (s) ,  one  can  solve  the 
nonhomogeneous  problem: 

i  =  P(t)l^t(t),  iUo)  =  lo,  = 

Its  solution  has  the  form 

I 

<0 

There  occurs  the  following  dependence  between  the  fundamental 
solutions : 

I 

V{t„t)  =  \U{l,  t)df, 

ft 

or 


U{t„  =  (6) 

For  the  proof  of  these  relations  it  is  sufficient  that  the  expres¬ 
sion  in  the  right  member  of  the  first  equation  satisfies  the  second 
group  of  conditions  (4). 

a.  TFIE  REMAINDER  TERM  OF  COWELL’S  METHOD 

Let  us  consider  the  vector  function 

G{z)  =  {C?')(z),  G<*)(z) . Gt->(z)} 


of  the  scalar  argument  z  and  introduce  for  it  the  following  differ¬ 
ences  ; 


1 


G(1)-G(0).  Ci*>=C?-C<V 

t  t  B 

(;=1,  3, 


5  .  .  .). 


(7) 


0 


We  state  the  interpolation  formula  of  Sterling 
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Then  successive  integration  of  System  (lO)  gives 


(11) 


— Ai*=  J  dt'\  J  (/*+,  — O /(/')«//' =A*J(1  — «)/(/»-+- Ar)rf*. 

*k'k 


Writing  an  analogous  relation  for  the  points  ^t-i,  tk  and  com¬ 
bining  this  with  Equation  (ll),  we  obtain  finally 

1 


A’|(l  —  z)G(z)dz, 


(12) 


where  it  is  assumed 


^  (*)  =/ Uk  hz)-+-f  (ft  —  Az). 


(13) 


We  introduce  now  for  the  function  AO  the  differences,  analogous 
to  (7): 

/k=/(tk), 

Cj.  =/•«-/• 

^  * 

A'=C±-A'L± 


Then,  noticing  that,  according  to  (13), 

c,'*-'’  =  2/i’^’,  =  0  (/  =  1,  2,  . . .). 

and  substituting  (s)  in  (12),  we  obtain 


(14) 
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where 


tty  =  2 j  (1  —  (7  =  1,  (l5) 

0 


If  the  function  m  has  a  continuous  derivative  of  order  2X-+-3, 

1,  according  to  (9),  the  remainder  term  can  be  presented  in  the 
form 


^  _  1 
Vt  —  2 


1  _  rf**+3(5(C,) 

2  “ti+s  3J5ITS  * 


(16) 


To  be  convinced  of  this,  we  introduce  the  numbers 
/n„  A/„  /  =  !,  2,  ...  n,  representing  respectively  the  minimum  and 
maximum  (2v-4-3)th  derivative  of  the  function  G^^(z)  in  the  intei*val 
[ — •*  —  1,  i'-4-  11: 


Further,  noticing  that  for  the  factor  (1  —  2)i4|^s(«:) 

represents  a  quantity  of  constant  sign  (we  say,  positive),  we  multi¬ 
ply  it  by  the  inequality  and  integrate 


[  (1  —  *) 


rfS'.+SCT'XCiCf 


Hence,  by  the  known  property  of  continuous  functions,  we  prove 
the  existence  of  a  point  C  in  the  interval  [ — > — 1,  >-*-1],  such 
that 


f . . 

J  (1  -  *)/l^+j(*) - - •  * 


1 

2 


— • 
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Integration  of  the  first  temi  from  (9)  is  carried  out  analogously. 
For  sufficiently  small  h  expression  (16)  can  be  transformed  to 
the  form 


Cti 


</*-**/ Ut) 


(16’) 


Equation  (14)  presents  Cowell’s  method  in  the  difference  form. 
Corresponding  to  it  is  the  following  sum  method 

j:  «■.] .  (17) 


where 


We  cite  a  few  values  of  the  coefficients  (l5) : 


1 

7 

12 

»s 

~  180 

1 

37 

~  240 

“6 

S040 

31 

199 

60480 

»7 

“  129600 

289 

1 

3628800 

3024 

3.  R0UNI>-0FF  ERRORS  AND  THE  ERROR  IN  COWELL’S  METHOD 
If  in  the  precise  formulas  (l4),  (l7)  the  remainder  terms 
are  omitted,  if  for  all  Xu  their  approximate  values  are  substituted, 
and  if  the  necessary  round-off  is  effected,  then  we  obtain  Cowell’s 
difference  method  and  sum  method,  respectively. 

We  shall  assume  that  System  (lO)  with  the  initial  conditions 

X{t^)—xt,  =  (18) 


12 


is  integrated  by  the  sSm  method  with  the  step  A '(within  the  domain 
where  the  sought  for  solution  exists),  and  as  a  result  of  the  integra¬ 
tion  we  obtain  the  table  of  values 


.  • Xf),  X\t 


X,, 


(19) 


at  the  corresponding  times 

•  •  •»  n  ^o>  ••••  •  •  •• 


SO  that 

We  introduce  the  quantity  t*,  defined  by  the  equation 

=  i  (^  =  0,  1,  . . .),  (20) 

/■J 

where  It  is  assumed 

$(/,  x)  =  h*F(t,  X),  =  'r)l’>  =  A*a>t,  etc. 

The  quantity  is  called  the  error  of  the  sum  method. 

If  System  (lO)  is  integrated  by  the  difference  method,  then  the 
quantity  t*,  which  is  defined  by  Equation  (20),  is  colled  the  error 
of  the  difference  method. 

The  error  arises  mainly  due  to  round-off  errors  at  each  step  of 
the  numerical  integration. 

For  investigation  of  the  question  of  calculation  of  the  error  we 
limit  ourselves  to  a  very  simple  case  of  Formula  (l?) : 

=  (21) 

If  one  digresses  from  round-off  errors,  then  the  computations  by 
the  method,  corresponding  to  the  precise  equation  (2l) ,  can  be  pre¬ 
sented  in  the  form  of  the  following  algorithm. 
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1.  The  approximation  of  the  values  of  the  initial  conditions  (l8), 
Xf,  X~i  are  given  and  we  enter  them  in  the  table  of  values  (19) ;  we 
compute  and  also  the  quantities 


2.  From  the  knovm  ^_i  we  extrapolate  (we  predict)  the  value 

of  <1),  we  compute  by  the  formula 


and  we  assume 


3.  We  compute 


and  we  assume 


4.  We  compute 


(p-i) 


;{  =(5(-*)  4._L  (P 

■^1,0  -  1  12  1. 0. 


I  —  1  12  1. 1 


and  we  form  the  difference 


if  it  is  found  that 
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where  X  is  some  preassigned  number,  then  we  finish  this  iteration 
process  at  the  first  step  of  integration  and  enter  the  value  of 
in  the  table  of  values  (19)  in  place  of  If,  however,  the  above 

inequality  is  not  satisfied,  then,  repeating  the  operations  analo¬ 
gous  to  those  set  forth  in  3,  we  construct  the  following  approxima¬ 
tions  : 

^l.t*  ^I.  S»  •  •  M  ^1,  I 

until  the  inequality  is  satisfied 

then  we  enter  the  vector  in  the  table  of  values  (19)  in  place 

of 

5.  If  in  the  table  of  values  (l9)  is  already  written  all  Xj  up  to 
7==^  —  1  inclusive,  then  for  computing  Xt  we  find  the  value  of 
defined  by  the  equations 

2* (82) 
~T  /■=» 


or 

02’) 

we  extrapolste  the  value  of  ,,  and  construct  the  following 
approximations  by  the  formulas 

=  /  =  o,  1,  ... 

The  iterations  are  continued  until  the  inequality 

— i+i  II  =  II  II  ^  1. 


(23) 


(33') 


is  satisfied. 
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Thus,  if  all  X*  from  (lO)  aro  obtained  in  /  approximations  on 
the  ^ th  step,  that  is  if 

X,  =  X,,, 


then  it  is  necessary  for  the  condition  (23')  to  be  fulfilled. 
The  quantity  is  called  the  error  of  method  (’2l) . 

[jOt  us  suppose  now,  that  all  ‘P*, /,  4'^  are  computed  with  the 
errors  /,  ft,  subject  to  the  inequalities 


(^>1,  />1), 

(^>-1), 


wliere  (i  is  a  fixed  number.  This  f’;enerates  in  the  quantities 
the  errors  (s^,  which,  accordin/^  to  (.22)  suid  (22'),  satisfy  the  rela¬ 
tions 


(24) 


Lot  us  suppose  further,  that  in  the  computation  of  Xii,i, 
round-off  errors  jJi*. i,  pi»,  a"<'  introduced,  satisfyinf'  the  inequalities 

Iliv.Kh 

II  f‘it  II  <  h  (k  >  1), 

whore  p  is  a  fixed  number  (for  example,  if  all  Xn.i,  X^  are  computed  to 
s  decimni  places,  th'^n  it  ff'llovjs  that  [i  —  ylO"'). 

With  regard  for  the  errors  introduced  above,  Equation  (23)  and  the 
condition  (23')  must  bo  written  in  the  form 

=  P»  '  or.  HM 

^ who re 

or.  error  in 

li'J’lr.  I  —  'I’lr,  l-H  •  fir,  I — fir,  (+1  ||  ~  ||  .^  1  (fe  ^  1). 

ir, 


(25) 


(25') 


Proceeding  now  to  the  calculation  of  the  error 

we  introduce  the  auxiliary  quantity 

t; -Jf, s'f’*  (•) 

then,  according  to  (25)  and  (25' ) >  we  find 

tj,  =  Pi  e  r .  ^  fit  n  P*'  “*■  ®  (n) 

(^  ^  !)• 

Hence,  observing  that 

t»=A*ti,  .  •(**) 

and  applying  the  operator  A*  to  both  members  of  the  above  equation, 
with  the  aid  of  (24)  we  find  finally 

where  it  is  assumed 

?»  —  ?»,  1+1. 

Remark.  If  one  makes  use  of  the  defJ.nition  of  and  intro¬ 

duces  also 

t 

then  Equation  (*)  for  *  =  0,  —1  gives 

Hence,  according  to  (**)  it  is  easy  to  show  that  Formula  (26)  can 
be  retained  for*»=o.  1,  if  Pji  P— i  eie  replaced  by  the  quantities 
Poi  P-ii  and  likewise  it  is  assumed  tnat 

a«  ■=  p_|  =  0,  Ifl  =  ^-1  =  h 

and  in  the  expression  er.  (]^)  considered  that 

=  t_|  =  0. 
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If  the  integration  of  System  (lO)  is  carried  out  by  the  difference 
method  corresponding  to  the  exact  formula 

then,  repeating  all  of  the  above  discussion  with  the  corresponding 
changes,  we  come  to  the  following  expression  fo”  the  error 

R  (n)  n  ^  (’26  ’ ) 

where  by  p,j^  are  meant  the  round-off  errors  for  subject  to  the 

condition 


\\K\\<v-'  (^>0). 

and  all  remaining  quantities  have  the  previous  meaning. 

Later  we  will  make  use  of  the  expressions  ('26)  and  ('26')  in  place 

of  the  approximate  values  for  the  error  defined  by  Equation  (20). 

Remark.  The  quantities  were  defined  as  errors  arising  from  the 
computation  of  the  vectors  if  in  the  latter  the  components  of  the 
vectors  Xi  are  considered  as  precise  numbers.  Following  Newcomb  and 
Brouwer,  we  will  consider  these  errors  as  independent  round-off  errors. 

4.  THE  LINEARIZED  ERROR  ESTIMATE  OF  COWELL'S  METHOD 

Let  us  suppose  that  the  times  ....  t-tt  •••«  ••• 

correspond  to  the  quantities  (l9)  obtained  by  numerical  integration 
of  System  (lO)  with  the  initial  conditions  (18)  by  Cowell's  method. 

This  means  that  these  quantities  satisfy  the  following  relations: 

=  %  {k>0),  (27) 


where  it  is  assumed 


Fk^F(lk,Xk),  =  etc. 
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and  by  ia  meant  the  error  which  is  computed  by  Formula  (26)  or 
('26’). 

Designating  by  JE(/)  the  exact  solution  of  System  (lO)  with  the 
initial  conditions  (18)  and  assuming  we  introduce  the 

following  definitions: 

a)  the  error  of  approximation  of  the  solution  Xt  (or  the  error  of 
Cowell's  method) 

=  Si*),  si"));  (28) 


b)  the  Jacobian  J{t,  Ji)  of  System  (lO) 

.\i  _  ^ti),  ^(») . j<«) 

(here  jE  represents  the  set  of  the  independent  coordinates) ; 

c)  the  Jacobigin  in  the  unknown  exact  solution 


y(o-y(^  i(0); 


d)  finally,  we  obtain 

1 

Qt  =  Jy(E»,  ^»-4-aS»)</a 
0 

[The  quantities  Qt  defined  only  for  the  condition  that  the  segments 
Lt,  connecting  points  (/»,  ^»)and  (^.  ■.^»)  of  theCn-*-!)  dimensional 
space,  are  contained  in  the  region  of  definition  of  the  solution 

.e  (0  ] . 

Subtracting  (l4)  from  (27)  and  using  the  obvious  identity 

Ft — /»=QA. 


we  find 

* 

y-i 


(29) 
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where  is  the  central  difference  operator  of  order  2j. 

Equation  (29)  is  called  the  equation  of  errors. 

In  order  to  simplify  the  problem,  we  replace  this  equation  by 
the  following  approximation; 

^  —  —  (29* ) 

Equation  (29’)  is  called  the  linearized  equation  of  errors 
(Lozinsky) . 

Let  us  introduce  the  fundamental  solutions  mu,  0,  y(t», Oof 
the  system  (f=y(0ff,  defined  by  the  conditions  of  type  (4) .  Then 
the  general  solution  of  Equation  (29')  by  analogy  with  (5)  can  be 
written  approximately  in  the  form 

0-^- (^>1).  (30) 

where  it  is  assumed 

(31) 

(32) 

(33) 

Every  estimate  according  to  the  absolute  value  of  the  right 
member  of  Equation  (30)  is  called  the  linearized  error  estimate  of 
method  (14)  or  (17). 

The  quantity  ircDi’Csents  the  error  created  by  the  errors 

1 — 1 

in  the  initial  conditions  (I8)  and  is  called  the  error  of  the  ini¬ 
tial  displacement. 
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Kk  =  U{t„  to  \  -  8_,), 

t-i 


Ir-l 


=  —  ^V{tm,  tO^rnh. 

msO 


»  I 


The  quantity  S*,  t  in  the  renult  oT  round-off  at  each  step  (and 
also  is  due  to  the  error  of  the  integration  method)  and  is  called  the 
round-off  error. 

Finally,  the  quantity  is  the  result  of  omitting  the  remainder 

term  in  the  integration  formulas  (14)  and  ('17)  and  is  called  the  quad¬ 
rature  error. 

Substituting  (16*)  in  (33)  and  approximating  the  sum  by  an  inte¬ 
gral,  we  find  the  estimate  for  the  quadrature  error 

1 ,  I  <  f  V{^,  t,)  dl  + 

>0  (34) 

^i» 

or  disregarding,  for  sufficiently  small  h  ,  the  second  term  of  (16') 

K jV(E,  f,)  ■  dL  (34 ' ) 

Usually  in  practice  the  quantities  T<,hs{%)6)  are  limited 'hy- 
the  absolute  value 

Taking  absolute  values  of  (32)  and  replacing  all  by  t,  we 

obtain  the  estimate  of  the  round-off  error  in  the  form  of  Lozinsky 

(35) 

where 

2  1^0 (5.  Ml <^5. 

i. 
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According  to  (26)  and  (26’).  the  sum  method  t  takes  the  value 
T  — •jp-»-4  er. 

and  in  the  case  of  the  difference  method,  the  value 

>s-ip  +  4er.(i)M+|.’H-ix, 


where  it  is  assumed 

Estimate  (35)  can  be  improved  a  little  if  use  is  made  of  the 
following  identity 

2  V{1„,  =  t,)-2V(t„  M]6«-+- 

inbO 

(*) 

rnml 

where  8^  are  arbitrary  vectors.  Introducing  the  definition 


and  utilizing  the  relations 

tu)-2V{u,  /»), 

V"(/»,  /»)  =  0,  h)f^hl, 

easily  deduced  from  (3),  (4),  and  (6),  we  shall  rewrite  the  identity 

(*)  in  the  form 
»»1 


2  2V(U,  Mi 6,  -4- 


/•*2  h)K-^hK 

fdsl 
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We  shall  nov;  substitute  (26)  in  (32)  and  the  result  of  the  sub¬ 
stitution  we  shall  write  in  the  form 

Sjt,  t  =  Sv.  p  -4-  p.  -»  p  -4-  X,  (36) 

where  it  is  assumed 

Sm=  S  V{t.,  tk)T, 

ms=0 

»->  n  (37) 

V{t„, 

m=0 

Let  us  assume  in  the  identity  (**)  6„  =  (i„  aiid  make  with  its  help 

the  estimate  replacing  all  by  their  maximum  value  (*; 

then  we  shall  obtain 

whore 

u 

The  estimate  of  the  temis  p.,  S^.x  from  (36)  is  made 
analogously;  thus  it  is  revealed  that  (at  least  for  equations  of 
celestial  mechanics)  all  these  terms,  differing  from  the  term 
by  the  small  factor  A’,  represent  small  values  in  comparison  with 
it,  so  that  the  estimate  (35)  for  the  sum  method  can  be  written  in 
the  form 
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(38) 


As  for  the  estimate  of  round-off  error  in  Cowell’s  difference 
method,  substituting  (S6')  in  (32)  and  reasoning  as  in  the  preceding 
case,  one  can  show  that  here  the  principal  terms  will  have  the  values 
P  and 

(37’) 

™=(»  n 

and  with  respect  to  them,  the  terms  p.,  S»,x  will  give  small 
corrections,  so  that  the  estimate  (35)  can  be  replaced  by  the 
following: 

1  Sf,’.  I  I V,  +  ‘l'.V  I  <  M" .  (38  • ) 

Thus,  knowing  the  fundamental  solutions  0*  ^(^o»  0  of  the 

system  P  —  J{t)S  and  using  the  estimates  (34')  and  (38)  or  (34')  and 
(38'),  we  can,  according  to  (30) ,  estimate  the  error  (28)  in  the 
integration  of  System  (lO)  by  Method  (17)  or  (14) : 

1 8*  I  ^  I  St,  0 1  -+- 1  S»,  j  I  1 5)r,  p  I , 

Remark.  By  integrating  the  equations  of  celestial  mechanics 
always  with  a  sufficiently  large  number  of  differences  in  the  integra¬ 
tion  formulas  and  with  sufficient  accuracy  in  the  initial  conditions, 
so  that  in  the  right  member  of  the  above  inequality  the  first  two 
terms  have  practically  no  influence  on  the  estimate  of  the  error 
then  this  inequality,  taking  account  of  (38)  and  (38'),  can  be 
written  in  the  form  ,  ... 

(39) 

(for  th.,  BUT,  method)  or 

(for  the  difference  method)  where  the  quantities  are  taken  from 

(35). 

The  estimates  (39)  and  (39')  can  be  improved  if,  as  adopted  in 
celestial  mechanics,  round-off  errors  are  considered  as 

independent  random  variables  in  the  sense  of  the  theory  of  proba¬ 
bility  (see  the  remark  at  the  end  of  Paragraph  3). 
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5.  ON  CONDITIONS  OF  'ITIK  APPLICABILITY  OF  HIE  NORMAL  DISTRIBUTION 
LAW  OF  PROBABILITIES  TO  THE  ESTIMATION  OF  ROUND-OFF  ERROR 

Some  real  variable  E  defined  by  a  random  quantity  is  called  sto¬ 
chastic  (or  simply  "random  variable"),  if  for  any  fixed  x  the  proba¬ 
bility  of  the  inequality  5<[x,  designated  by  <p(x),  is  known;  by  varying 
X  from — cn  to  -i-oowe  obtain  the  function  'p(x),  which  is  called  the 
distribution  law  or  the  distribution  function  of  the  random  variable 

The  relation  of  the  random  variable  {  to  its  distribution  function 
<p(x)  is  written  in  the  following  form: 

P(l<x)  =  cp(x). 

The  argument  of  P  can  be  transformed  to  equivalent  inequalities 
without  changing  their  probability. 

Axioms  of  probability  impose  on  the  distribution  function  of  any 
random  variable  certain  limitations: 

1)  <f{x)  monotonically  increases  in  the  interval  — oo  < x ! 

2)  <p( — oo)  =  0,  <(>(-+- oo)=l,  etc. 

Knowing  the  distribution  function  of  the  random  variable  5,  for 
the  given  quantities  x,,  x,  one  can  find  the  probability  of  the 
inequality  x,  <C5<Cxj  by  the  formula 

P(x,  <  5  <  xj)  =  <p  (xj)  —  f  (x,), 

whence,  in  particular,  it  follows  that 

P(j^f<x)  =  f(x)  —  f(—x).  (40) 
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.F 


The  random  variable  C,  having  the  distribution  function 

9  1  . 

is  called  normal,  and  the  function  C(z) 

—  » 

called  the  normal  distribution  law  of  probabilities.  According  to 
(40),  for  the  normal  random  variable  C  we  have 

P{|C!<z)=<I.(z)=-j|rJr'^'‘rfz.  (41) 


We  cite  a  few  values  of  the  function 


t 

fU) 

z 

•M*) 

0.03 

00239 

0.674 

0.5000 

0.1 

00798 

1.0 

0.6826 

(42) 

0.2 

0.1586 

1.5 

0.8650 

0.3 

02360 

2.0 

0.9544 

0.4 

03108 

3.0 

09973 

0.5 

03830 

4.0 

• 

0.9999 

For  the  random  variable  i 

,  the 

distribution  function  of  which 

has 

a  continuous  derivative  p{x), 

the 

following 

definitions  are  intrO' 

- 

duced: 

1)  expected  value 


■♦-ao 

A'  5=  j  xp{x)dx^a\ 

-a>  • 


2)  variance 


D{\)~  EiS-af^b; 

3)  third  absolute  central  moment 

/r|5— a|’  =  c. 

If  1  Is  a  constant  quantity,  then 


=  D{\,  5)=X*£>(I). 
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The  random  variables  are  called  totally  independ¬ 
ent  if  the  probability  of  the  inequality  for  any  fixed  mis 

not  dependent  on  the  values  the  remaining  random  variables  take. 

THE  CENTR/VL  LIMIT  THEORM  OF  THE  THEORY  OF  PROBABILITY 

Let  5*  be  some  totally  independent  random  variables, 

possessing  the  expected  values  E  =  the  variances  b„  and  the 

third  absolute  central  moments  c„(l^/n^A:). 

Then  the  distribution  law  of  the  quantity 

k 

r  _  m=8l 

k 

where  Bk^^b„  tends  to  the  noiroal  law  G(r)  for  ^-*00,  if  is  satis- 

m=l 


fied  the  following  condition  of  Lyapunov; 


«|»»  =  ^^^-5 - ►  0  for  k-*  CO. 

b1 

The  theorem  signifies,  that 


C(x), 


or,  according  to  (4l) , 


Returning  now  to  the  notation  of  Paragraph  4,  we  ahall  consider 
the  round-off  errors  from  (37)  as  independent  random  variables 

having  the  same  distribution  functions 

ip(x)  =  0  for  x<  — p, 

'f  W- for-  —  P “^P. 

<p(x)  =  1  for  >  P 

(uniform  distribution  of  a  random  variable). 

The  quantity  p,  appearing  in  the  definition  <p(x),  must  satisfy 
the  inequalities 

(43) 


Applying  the  theorem  to  the  random  variable  we  find  succes- 

si  Vfily 


4-  .2 


-r^O  /=l 


where  it  is  assumed 


(44) 


I,  /’> 


The  above  theorem  leads  to  the  following  statement.  If  the  con- 


WR  opply  the  normal 


in  fulfilled,  then  to  the  round-off  error 
distribution  law  of  probabilities 


(45) 


Remark.  In  all  cases,  which  we  subsequently  encounter,  the  con¬ 
dition  (*)  is  fulfilled,  and  we  will  not  dwell  on  it  further. 


Let  us 


introduce  the  coefficient  of  overestimating  of  Formula 
which  is  itself  a  random  variable 

P(v<!/)=p(j  VW <  I  SI?, !)  =  1  -  (f )  . 


(45) 


The  quantity  </,  defined  by  the  equation  rf~P(Yj<l)  =  l  — 
represents  the  probability  of  underestimating  and  must  be  considered 
as  the  defect  of  Formula  (45)  for  a  given  z. 

Of  greatest  interest  is  the  quantity 

p(!/)-  P(1<t,<!/)  =  .1.  (z)  -  «1)  (-i)  ,  (46) 


which  characterizes  the  probability  of  "cleanly"  overestimating  given 
by  Formula  (45) . 

For  z=:3  Formula  (45)  gives 

with  a  probability  of  0.9973,  where  the  quantities  P,  are  deter¬ 

mined  according  to  (43)  and  (44). 

With  regard  for  the  observation  made  at  the  end  of  Paragraph  4, 
one  can  write  the  following  probability  estimate  for  the  error  of  the 
sum  method  (l7) : 


The  probability  of  underestimating  will  be  t/ =  0.0027. 

With  the  help  of  (42)  and  (46)  one  can  determine  the  limits  of 
"cleanly"  overestimating,  given  by  Formula  (47)  or  (48)  and  the 
probabilities  of  their  realization. 


y 

ply) 

6 

0.6143 

10 

0.7613 

30 

0.9175 

100 

0.9734 

Thus,  for  example,  according  to  the  resulting  data,  one  can 


assert  that  with  a  probability  of  0.7613  Formula  (47)  gives  a  "clean" 
overestimate  by  less  than  a  factor  of  10. 


Remark.  Repeating  the  above  discussion,  we  can  write  the  proba¬ 
bility  estimate  for  the  terms  from  (36)  if  the  errors 

are  considered  independent  random  variables  and  use  is  made 
of  the  identity  (*)  from  Paragraph  4  [here  are  obtained  small  correc¬ 
tions  to  the  right  members  of  the  inequalities  (47)],  and  also  the 
estimate  of  the  tern  (37')  if  the  errors  are  considered  inde¬ 
pendent  quantities.  Hence,  making  the  estimate  of  (-37')  together  with 
(37),  we  shall  obtain  the  approximate  estimate  of  round-off  error  for 
the  difference  method  (l4)  _ 


v^3 


£ 


which,  according  to  the  observation  at  the  end  of  Paragraph  4,  can  be 
extended  then  to  the  error 


6.  ESTIK\TI0N  OF  THE  ERROR  OF  NUMERICAL  INTEGRATION  OF  THE  PROBLEM 
OF  TWO  BODIES  (ELLIPTICAL  MOTION) 


Let  us  consider  the  system  of  equations  of  undisturbed  motion 


,.3  — ^  ,.1  — 0* 


(49) 


r*  =  X*  y*  +  r* 


and  let  us  assume  that  it  is  integrated  by  the  sum  method.  The  gen¬ 
eral  solution  of  System  (49)  has  the  form 

X  =  a  [P,  (cos  £  —  e)  v/ 1  —  e*Q,  sin  £], 

y  —  a  [P,  (cos  £  —  e)  v/l  —  e*Q,  sin  £],  (50) 

z  =  a  [P,(cos  £  —  e)-¥  —  e*Q,  sin  £], 
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where 


E  —  e sin E  =  n(t  —  7^,  n  =  Ka  *  ’ 


Here  P^,  P^,  P„  Q,  denote  direction  cosines  (Subbotin, 


1941) 


Let  us  introduce  the  vectors 


x  =  {x,  y,  z)^{xi'\ 

x*{t)^{a(cos  E  —  c),  flVl  —  c*  $in£,  O) 


and  the  orthogonal  matrix 


S--\\s,j  11  = 


P.,  Q„ 


P„  0„ 

P„  Q..  ^\/l-Pi-Qj 


(50') 


(51) 


The  signs  in  the  third  column  of  the  matrix  S  are  taken  in  such  a 
way  that  the  orthogonality  condition  is  fulfilled,  i.e. , 

SS  =  L  (52) 


Then  Formula  (50)  can  be  written  in  the  form 

je(/)  =  5je*(0. 


In  the  notation  of  Paragraph  4  we  find 


Jit.  =  ^ 


3xz 


3x:*—  r\  3xy, 

3->ry.  3y*— P,  3yz 


3jrr, 


3yz,  3z’  —  p 


(53) 


Desiring  to  simplify  the  algorithm  for  computing  the  fundamental 
solutions  U{t„,  t),  V{t„  t)  of  the  system 

S~J{t)y, 
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we  introduce  the  matrix 


J*{t)  =  J{t,  x*{t)) 

and  the  fundamental  solutions  t),  0  of  fhe  system 

r  =  /*(0.9*. 

Then  according  to  (l),  (2),  (52),  and  (53),  we  find  successively 

(/({„  f)  =  S[/*(f<„  0^, 

V(t„  ()^SV*(i„  f)S. 


Proof  of  the  relations  (54)  is  immediately  and  sufficiently  con¬ 
vincing  by  the  fact  that,  using  (53),  their  right  members  satisfy  the 
conditions  of  Type  (4). 

Matrix  J*{t)  has  the  form 


/•(/) 


A,  B,  0 

B,  C,  0  , 

0,  0,  D 


whe  re 


Hence 


*  „  j(3  —  e*)  cos*  E  —  4e  cos  E  ■+■  (3e*  —  1) 
'  ”  (1  —  e  cos  E)^ 


u _ j3 VI  —  e*  (slnfi'cos^’  —  csInE) 

(T^cosE)»  > 

^  _  i  (2f*  —  3)  cos*  E  2e  cos  E  {2  —  3«*) 
~  ^  (1  —  e  cos  £)' 


D  =  a» 


-1 

(1  —  cos  Ey  * 


<5(^.  t)  ^0, 
0"  0,  0="0, 
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and  the  quantities  v*,,  u,  ,,  v[,,  wj  j  satisfy  the  equations 

:  /!«;  ^  -+  Bv]  ^ 

(y  =  l,  2). 

Dvl, 

The  solution  of  these  equations  has  the  form 


< .  0  =  „(1-.cos£H1-«co7£;;)  sin  £  cos  £  C  sin  £|  X 

X  cos’*  £o  [(c*  —  1)  cos*  £  «-  (e’  —  e)  cos  £  -»-  2  (1  —  c*)]  cos  £9  sin  £9  -+- 
-4-((e  —  c®)sin£cos£-t-2(e*-<-  l)sin  £]cos  £9-»-[ — ccos*£  — 

—  2(c*-4-  l)cos£  K  5e)sin  £9  +  1—2(1  —  c*)  sin  £cos  £  — 5e  sin  £]  — 
—  3  sin  £  sin  E^(E —  £9)): 

( |cos*  £  ♦-  c  COS  £  -  2J  cos*  £9  -4- 

(sin  £  cos  £  c  sin  £)  sin  £„  cos  £9  -4-  ( — e  cos*  £  +-  2  cos  £  —  c]  cos  £9  -4- 
-4-( — csin£cos£  I  "2  sin  £]  sin  £9  4- [cos*  £  4  ecos£ — 2]  1 

•  3sin£cos£9(£  —  £9)): 


v/ 1 _ e- 

t'j.  1  (^1  0  —  „  (1  —  e  cos  /■)  (I  —  e  cos  £9) 


-COS*  £ -4- e  COS  £ — 1]cos*£9  t- 


-t  [ — sin  £  cos  £  -4  c  sin  £]  sin  £9  cos  £9  4^  [ — c  cos*  £  — 2  cos  £  —  c]  cos  £9  -f- 
4-[ — esin£cos£  —  2sin  £]  sin  £9 -4  [2cos*£  4  ccos  £-4-2]  4 
4-  3cos  £  sin  Eo(E —  £9)); 

<  0  =  ^-7c7s~g)(T-;  c75£9)  ^  ^  -  e  sin  cos*  £9  -4- 

4-  [ —  COS*  £  -4-(c*  4-  e)  cos  £  —  I]  sin  £9  cos  £«  4  [ —  (c*  4  e)  sin  £  cos  £  4  -  2  sin  £]cos  £9  -+- 
I  [ccos*  £  —  2  cos  £  4  cj  sin  £9  4  [sin  £  cos  £  —  c  sin  £]  —  3  (1  —  e*)  cos  £  cos  £9  (£  —  £9)); 


wj  j(/9,0  =  —  [sin  £  cos  £9  •  (-  cos£-(  c)sin£9 — csin£]. 


The  correctneaa  of  Equations  (55)  is  ascertained  by  direct  sub¬ 
stitution  of  their  rifijht  members  in  the  above  equations  and  verifica¬ 
tion  of  the  appropriate  initial  conditions  [see  (4)].  Here  by  £9  is 
meant  the  initial  value  of  the  anomaly 

£9  — esin  £„=  T). 
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PassinR  now  to  the  definition  of  the  matrix  U*  /),  we  have 
according  to  (6) 

( 1  _  e  cos  £■)  (I  —  e  cos  foP  I  ^  ® ^  —  c)  (£  —  Eq)  -I-  .  .  .  ), 

<*(<«-  ^)=^(l-ocos£m-^cos£oy  !3  Sin  £  sin 

0 •■=  (Tzr,-7os £) 0 - « cos goP  ^ £o - c) (£-£,)  +  ■■  ■ ), 

5  ('« .  0  =  (i-cVosgrii  -  I -3 (1  -  «')  sin  £o  cos  £  (£  -  £o)  -4-  . . , 

^n.  3  (^ii>  O  '"  •  •  •  I 

where  the  dots  denote  periodic  terms,  which  for  us  have  no  importance. 
The  relations  (54)  give 

'Sy.  2  *~^3, 3(^n>  0  S(,  3  Sy,  3 

and  the  analogous  equation  for  0> 'thence,  introducing  the 

quantities 

sin  £^—  \/l  —  c*s,,2Cos£», 

,  (cos  £»  —  e)-i-\/l  —  3  sin  £*  (/=!,  2,  3), 


where 


£»  -  e  sin  /;V  «  (^  — •  T"), 


according  to  (55)  and  (56)  we  find  to  within  the  secular  terms 

«’<.  y  {in<  ft)  =  ;r(r-77^wi7y (T^t^osIQ 

'‘>=-rrrTco;7:;nT^  ^o). 


(57) 


Substituting  (57)  in  (3I)  we  obtain  the  expression  for  the  error 
of  the  initial  displacement 

.  ,'X> .  rX'hti  -  . 
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Equation  (34')  with  regard  for  the  relations  (53)  and  (54)  gives 
the  estimate  for  the  quadrature  error 

I, 

Hence  with  the  aid  of  (50*)  and  (55)  we  find  to  within  the  secular 
terms  (for  simplicity  the  fonnulaB  are  written  for  c  =  0) 

£o). 

We  turn  now  to  the  calculation  of  the  quantities  Nli*\  defined  by 
Equations  (44)  and  appearing  in  the  estimate  of  the  round-off  error 
(47).  For  this  we  notice,  that  these  quantities  fonn  the  main  diago¬ 
nal  of  the  matrix  VV,  which,  according  to  (54),  can  be  represented  in 
the  form 

VV  =  SV*  \?*S. 
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Elementary  calculations  according  to  (55)  lead  to  the  following 
expressions  for  the  quantities  which  for  simplicity  were 

written  for  e  =  0): 

Ar;n.i)  =  ^  j 3  sin‘  £*  (£»  -  £«)’  +-  6  sin  £*  cos  £*  (£*  —  £,)*  -i- 

sin*  Et  12  sin*  Etcos  (Et  —  £o)  -H  12  sin  £*sin  £ftj(£t  —  £o) 

-♦-^24  cos*£t  sin(£»  —  £#)  —  32  sin  (£»  —  £o)-t-  3(sin£*cos  £*  —  sin  £3003  £j)  — 

—  j(sin*£t  -«-  l)sin2(£t  — £o)  J| , 

*1  =  -^  {— 3  sin  Et  cos  £»  <£»  —  £„)»  -  3  (cos*  £»  —  si n*  £*)  (£*  —  £3)* -H 

cos  £»  —  12sin£t  cos£tcos(£i;  —  £3)  —  6  (sin  Etcos  E^-*-  sin  £3  cos  £t)J  X 
X  (£t  —  £3)-+- 1^24  sin  Et  cos  £»  sin(£k  —  £3)  —  y  (cos*  Et  —  sin*  Et) 

-s-  y  (cos*  £3  —  sin*  £3)  -+■  y  sin  EtcoaEt  sin2  (£»  —  ■^o)J| » 
yv;(*'  ”’  =  -^  {  3  cos*  Et  {Et  -  £3)*  -  6  sin  £t  cos  Et  (Et  —  E^f  -h 
j^y  y  cos*  Et  *♦- 12  cos*  £t  cos  {Et  —  £3)  12  cos  Etcos  £31  (£»  —  £3)-^ 

-*-j'24  sin*  Et  sin(£t  — £3) —  32  sin  {Et —  £3)  —  3(sin£t  cos  Et  —  sin  £3 cos  £3)  — 

—  y  (cos*£t-«  1)  sin2(£t  —  ^0)  j|  * 

^>=-^(j{Et-Eo)~j  sin  {Et  ~  £3)  cos  (£»  -  £3))  . 

Hence,  having  used  the  quantities 

<4'*  =  s.  ,  sin  Et  —  s.  ,  cos  Et 
:  *  (»  =  1,  2,  3), 

Y<*>  =  s,  ,  cos  Et  -+-  Sf  j  sin  £j 

from  (58)  we  find 

=  6a<;V;’(£»-£o)'  H[L3_6s?.3-»-y^'''Vl2o<;'*co8(£»-  £3)- 

-I-  12oi'W„‘)  1  {Et  -  £3)  -f-  [-8(1-  s* , )  sin  {Et  -  EJ  -  24^^)*  sin  (£^  -  £,)  -+- 
3  (o^*)4o_ ^  I  -  y  -H 1  s*  ,)  sin  2  (£* -  £3)])  (/  =  1,  2.  3). 


(59) 
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For  *7^0  ,  Equations  (59)  must  be  replaced  by  the  following 
(for  = 

'VP  =  «in  _ 9. ,in  (£»  -  £^* 

+  [4  -  •*>*i»’* 2^0  -♦*  (1’-  •*)  [24 (cos  £*  H-  .)  ay)*— 

- 12  sin  £,oy')YV)]  cos  £,  +  (1  -  **)  [6  (.*  4)  sin  Etof' 

-♦- 12  cos  £toy)ify)]  sin  £#  -4-  j|2«*  cos*  £*  5«  cos  £»  •+■  (7  —  y  **)!  7*“ 

9«*)  cos*  Ek  -t-  (4«  —  8e*)  cos  £»  -+  (  —  j  e*  —  T  **  "*“  ^4^  j 
[8e*  sin  £»  cos  £»  -♦•  (8e*  ■+-  4«)  sin  £»J  0^)11/)  j(£»  —  £b)  -•-  [ •  •  •  j|  1 

where 

s<,.i  sin  El  —  \J\  — e*s«,  1  cos  £i  1 
=  e*)cos  Ek  —  2e]s«,  1  -4-  \/\  —  e’ S(,  j  sin  £». 

Here  the  dots  denote  periodic  terms  differing  by  a  small  quantity  from 
corresponding  terns  of  (59). 

Inversion  of  Keppler's  equation  [see  (50)  ]»  gives 

E—  M-t-e  s\n  M-+-  4  c'  sin  2Af  +- . . 

where 


hence 

Ek  —  £#  =  f^fth  -•  2c  sin  cos  ^  4  c*  sin  knh  cos(M^  »-  Af^)  . 

I  f  corresponds  to  an  integral  number  of  revolutions,  then 

Ek —  E^  —  knh. 


(60) 
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I! 


For  a  sufficiently  large  value  of  knh  and  for  Formulas 

(47)  and  (60)  give  the  following  approximate  expression  for  the 
round-off  error: 


(61) 


where 


oV*=H  1  sin  Ek  —  v/1  —  c*  S(.  a  cos  Ek- 

One  can  compare  this  result  with  the  well  known  estimation  of 
Newcomb,  which  in  our  notation  can  be  written  in  the  form 

£ 

i|S».pii<0.225pA:*. 

From  the  standpoint  of  the  theory,  stated  in  Paragraph  5,  the  • 
probability  of  such  an  inequality  will  be  equal  to  0.1780  [see  (42) 
and  (46)  ] . 

If  use  is  made  of  the  estimate  of  round-off  error  in  the  form  of 
(35),  then  we  obtain  the  following  results: 

(62) 


Remark  1.  Let  us  consider  the  case  of  circular  motion  described 
by  the  equations 

=  0  (/  =  !,  2.  3). 

3 

2  —  const , 

i,  1 

where  n*  is  a  constant. 

Let  us  suppose  that  this  system  is  integrated  by  Cowell's 
method.  The  general  solution  of  the  system  has  the  form 

^0  =  r  cos  n/ -4- sin  n/)  (i  =  1.  2,  3). 

Calculating  the  estimate  of  the  error  of  numerical  integration  of 
this  system  by  Formulas  (47),  we  find  successively 

sin  n  (t  —  tn)  ,  .  _ 

Vi, ({to,  t)  - - 21 .  T.(,  j  (tn,  t)  =0,  I  yii ); 


[/-jP  . - = - 

^  ‘l**=  - 3 —  y  knh  —  —  sin  2knh 


{nh)* 


nh 


k  . 
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This  result  was  predicted  by  Brouwer. 

Remark  '2.  Let  us  consider  at  the  seme  time  the  equations  of  disturbed 
and  undisturbed  motions 

-f*  j 


‘R. 


K*jc 


=  0. 


where  R  is  the  disturbing  function.  We  shall  denote  respectively  by 
X„  ^the  results  of  integrating  these  systems  by  Cowell’s  method  with 
the  same  initial  conditions.  Then  because  of  the  small  value  of  the 
disturbing  function/?  one  can  put  —  XtvHl—  x  and  estimate  'the 
error  of  numerical  Integration  of  the  system  of  disturbed  motion  by 
the  formulas  of  Paragraph  6.  (This  hypothesis  is  a  supposition  of 
S.  M.  Lozinsky.) 
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TIIE  ACCUI.IULATION  OF  KUTffiRICAL  INTEGRATION  ERRORS  IN  SOME 
PROBLEMS  OF  CELESTIAL  MECHANICS 

A.  S.  Sochilina 

The  results  of  the  application  of  estimetions  of  errors  in  numeri¬ 
cal  integration  (Myachin,  1959)  to  numerical  examples  are  presented. 

Recently  in  connection  with  the  rapid  development  of  computing 
equipment,  methods  of  numerical  integration  for  the  solution  of  the 
problem  of  n  bodies  have  become  most  effective.  Such  enormous  works  , 

of  calculation  as  "Coordinates  of  the  Five  Outer  planets  1653  -  2060", 

"Coordinates  of  Four  Minor  Planets  1940  -  1960",  and  others,  were 
accomplished  by  methods  of  numerical  integration.  However,  the 
accumulation  of  errors  in  integration  essentially  reduces  the  quality 
of  the  numerical  methods. 

Errors  arise 'due  to  limited  precision  of  the  calculating  machine  I 

(round-off  errors),  neglected  differences  in  the  integration  formulas,  j 

and  inaccurate  values  of  the  initial  conditions. 

Errors  due  to  neglected  differences  can  be  reduced  to  a  minimum  by 
a  suitable  choice  of  the  interval  and  the  number  of  terms  in  the  inte-  i 

gration  formula.  The  calculation  of  errors  in  the  initial  conditions 
does  not  present  great  difficulty:  if  the  initial  conditions  are  com¬ 
puted  with  the  precision  with  which  the  computations  are  carried  out, 
then  one  can  consider  these  errors  as  round-off  errors  in  the  first 
step  of  integration. 

In  the  majority  of  cases  it  is  important  to  know  the  dependence 
between  the  number  of  steps  and  the  number  of  lost  digits,  in  order 
to  provide  for  the  necessary  precision  in  advance  or  determine  with 
what  error  any  or  all  quantities  were  obtained. 

With  this  aim  B.  F.  Myachin  (1959)  investigated  the  formulas  for 
estimation  of  the  accumulation  of  round-off  errors  in  numerical 
integration  of  the  equations  of  motion  of  the  two-body  problem. 
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If  the  eccentricity,  «,  ie  neglected  and  the  time  error  at  the  feth 
step  is  indicated  by  ,  we  can  write  these  formulas  in  the  following 
form: 

(k=\,  2,  3...;  /=1,  2,  3), 

where  ^  is  the  number  of  steps,  /is  the  number  of  the  coordinate 

(■*■.  y,  «). 

(nh)  '* 

/V(‘)  =3  3o(;)* (£,-£,)=>  4  £^)*  >  |^^_6s»3-4-yoi')*  4  12oi.0^  COS  (£,-£„)  4- 

12o(;)o(o]  {E,  -  £„)  4-  [-8(1-  si ,)  sin  (£,  -  £„)  -  24o(0»  sin  (£,  -  £„)  -4-  ( l) 

•4-  3  3  a<0=  _  3.  4-  1  s» ,)  s in  2  (£,  -£„)], 

<4‘>  =  s,_,  sin£^— s^^^cos  £,. , 

Ti'’  =  \iCos£^  4-s,  .,sin£^, 

S(,i.  S(,n  are  the  direction  cosines  (in  general  denoted  by 
P,  Q,  R),  n  is  the  diurnal  motion  of  a  body,  h  is  the  integration 
step,  £*  is  the  eccentric  anomaly  (in  the  given  case  the  mean,  since 
e==0),  and  p  is  the  maximum  round-off  error  in  the  computation  of 
the  right  members  of  the  equations  at  each  step. 

Vte  used  the  formulas  in  this  form  for  numerically  comparing  the 
error  8H),  obtained  by  numerical  integration,  with  the  predicted 

‘V>. 

1.  We  shall  carry  out  the  indicated  comparison  in  the  exsmples, 
which  were  specially  computed  for  this  purpose  (the  examples  were 
computed  on  the  electronic  machine  bESM  AN-SSSR) .  The  problem  of 
plane  undisturbed  motion  is  solved  viith  various  initial  conditions 
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(orbital  elements).  The  value  of  the  interval  of  integration  and 
the  elements  are  aeleeted  in  sueh  a  way,  that  one  revolution  is 
aeeomplished  in  exaetly  100  steps.  In  all,  three  examples,  eaeh 
having  1100  steps,  were  eomputed;  the  initial  conditions  were  defined 
by  the  following  elements: 

Ist  example  '2nd  example  3rd  example 

Mn .  0°  0“  0° 


Integration  step  A 


0°  0“ 

0°  0° 

0.04825380  0.04825380 

299'.' 128376  648" 

43!'35258794  20('0 


Integration  was  conducted  by  Cowell's  sum  method  taking  account 


of  fourth  differences 


^  240-^*^’  '  60480-^'^’’ 


where 


and  /Y  is  a  vector  with  the  components  X,  Y,  Z. 
The  results  of  the  integration 


Xu  X, . . 


were  compared  with  the  value; 


.^1,  .^2, 


computed  beforehand  by  the  formula  of  elliptic  motion  to  the  6th 
decimal  place,  and  the  difference 


A  —  Xt  =  Sn 


was  taken  as  the  net  accumulation  of  round-off  errors,  since  the 
influence  of  the  higher  differences  in  the  integration  formula,  in 
the  given  case  was  computed  with  the  maximum  precision 


In  the  future  we  will  not  take  these  errors  into  account,  since 
in  all  of  the  works  considered  their  influence  lies  beyond  the  limit 
of  precision  with  which  the  computations  are  performed. 

For  the  round-off  error,  f, ,  the  adopted  value  is  0.5  x  10“°;  that 
is,  the  precision  with  which  /  was  calculated  at  each  step.  One  can 
note  that  Formulas  (l)  are  considerably  simplified  if  the  estimates 
are  computed  for  the  points  ft  —  as  multiples  of  0,-^,  «,  ir. 

In  our  case,  since  the  problem  of  plane  motion  was  solved 
(s,,  =  S2i  =  ^»  Su  =  S2i  =  0)  integration  began  from  the  point  of 
perigee  (£o  =  0).  Formulas  (l)  are  simplified  still  further  and 
become 

5(01<40  =  -Jl^pv/yV(W  2;  ^  =  1,  2,  3  ...  ). 

{nhf 

/V<"  (£,)  =  £„  =  3£2  38E,  for  E,  =  2mir,  {m=  1,2. . .) 

/V'"(£»)  =  3£2  t  ]4£»  — 32,  for  £»  =  (2m  f-i)  k,  (2) 

/V<'>(£»)  =  ^  Et,  yV'*'(£»)  3El  -  lOE,  for  £»  =  (2m  -+- 1)  r, 

/V<'>(£»)  =  3£2  -«  14£»-»-32,  (E,)  =  j  E,  i- 8  for  £»  =  (2m-i 

Beginning  with  /fc>200,  the  ratio  is  of  the  order  of-i_: 

100  ’ 

therefore,  neglecting  the  first  power  of  £»  in  comparison  with  £J 
and  letting  En  —  nhk,  we  shall  obtain 


10.3pk^  , 
^^^')=3pk\ 


si*’  =  3pAr'‘  for  E It  — m-^, 

=70.3?Ar^  for  Ek=^{m  -i  y)’'. 


(3) 
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We  note  that  for  all  three  of  the  examples  under  consideration 
are  the  same.  This  is  stipulated  by  the  fact  that  the  product  nh  in 
all  the  examples  is  one  and  the  same  (w)  and,  in  addition,  the 
influence  of  eccentricity  is  excluded  in  the  cited  formulas. 

However,  the  eccentricity  does  not  very  strongly  distort  the 
result.  For  A: >200  the  formula  of  estimation  (Myachin,  1959)  is 
written  in  the  form 

l2 

*»  =  3p  I  Ot  I  -I - TT-  . 

In  Table  1  are  cited  the  computed  and  the  predicted  errors 

The  first  column  gives  the  mean  anomaly  for  all  three  "planets"*, 
the  '2nd,  the  appropriate  number  of  integration  steps;  the  3rd  to  8th 
columns,  the  true  errors  Sj*)  and  for  examples  I  to  III;  the  $th 

and  10th  columns,  the  values  of  and  the  estimate  of  the 

error  in  all  three  examples.  The  values  8^',  are  expressed  in 
units  of  the  sixth  decimal  place. 

Comparing  the  3rd,  5th  and  7th  columns  with  the  9th;  or  the  4th, 
6th,  and  0th,  with  the  10th;  we  can  judge  the  quality  of  the  obtained 
estimate.  Note  that  the  estimate  reflects  the  oscillatory  character 
of  the  accumulation  of  errors. 

'2.  From  the  standpoint  of  the  accumulation  of  errors  it  turned 
out  to  be  of  interest  to  exam.ine  the  coordinates  of  Uranus,  Saturn, 
and  Jupiter  obtained  by  D.  K.  Kulikov  by  integrating  the  Vlllth 
satellite  of  Jupiter  for  the  period  from  '24  January  1930  to  28  August 
1965.  The  integration  was  carried  out  on  the  electronic  calculating 
machine  BESM  with  steps  of  10  days  (all  1300  steps).  The  coordi¬ 
nates  of  the  plsnets  were  obtained  by  simultaneous  integration  of  the 
system  of  nine  equations;  the  initial  coordinstes  were  taken  from 
"Astronomicsl  Papers"  (l95l). 
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TABLF  1 


Since  the  cooi'dinateB  of  the  planets  published  in  "Astronomical 
Papers"  in  1951  were  computed  with  great  precision  (the  joint  system 
of  equations  of  motion  of  the  five  outer  planets  was  solved;  more¬ 
over,  the  calculstions  were  carried  out  with  14'  digits),  they  were 
taken  for  the  precise  solution  and  the  difference  jf*  —  Xk  was 
taken  for  the  true  error  of  the  estimated  coordinates.  The  compari¬ 
son  is  made  in  the  5th  decimal  place.  The  results  are  given  in 
Table  '2. 

For  an  estimate  of  this  error  by  the  formulas  of  B.  F.  Myachin, 
t  is  necessary  to  ascertain  the  error  made  in  a  single  step.  In 
the  given  case,  in  addition  to  the  round-off  error  p (p  —  0.5  •  10""), 
the  computation  of  the  right  members  of  the  equations  will  contain  an 
error  due  to  the  neglected  disturbances  from  Neptune  and  Pluto.  The 
magnitude  of  the  disturbances  amounts  to  1  •  10'";  that  is,  it 
exceeds  the  computing  error.  Of  the  three  planets,  Uranus  is  subject 
to  the  greatest  disturbing  influences;  since  during  the  investigated 
interval  of  time  Uranus  makes  only  a  half  revolution  (and  the  sum  of 
the  disturbances  from  Neptune  and  Pluto  has  a  constant  sign  during  a 
considerable  period  of  time),  while  Saturn  makes  1.5  revolutions  and 
Jupiter  3.5;  and  accordingly  the  disturbances  from  Neptune  and  Pluto 
nave  a  periodic  character.  For  the  estimate  of  the  error  due  to 
round-off  we  use  Formulas  (l)  (the  results  are  given  in  Table  '2  in 
columns  4',  7,  lO) .  For  the  estimate  of  the  error  due  to  neglected 
disturbances,  however,  it  is  impossible  to  use  these  formulas  since 
in  deducing  (l)  the  probability  law  of  the  distribution  of  random 
errors  was  used  snd  the  disturbances  do  not  obey  the  Isw  of  random 
errors.  The  only  acceptable  formula  in  this  case  is  the  following 
(Myachin,  1959) : 

(i  \,'2,3;k  1,2,3...),  (4) 
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TABLE  '2 


1  Number 

Round 

Error 

due 

Round 

Error 

due 

Round 

Error 

due 

Anomaly  I 

1  of 
steps 

XI,  — 

Off 

error 

to  neg¬ 
lected 
distur¬ 
bances 

Hit  —  Yy 

off 

error 

to  neg 
lected 
distur¬ 
bances 

n  -  A 

Off 

error 

to  neg¬ 
lected 
distur¬ 
bances 

Jupiter 


135 

90 

0 

0.1 

1 

0 

0.1 

0 

0.0 

ISO 

144 

0 

0.2 

1 

0 

0.3 

5 

0 

0.1 

2 

225 

198 

—2 

0.4 

0 

0.1 

0 

0.8 

210 

252 

—3 

0.5 

15 

-  1 

0.3 

3 

0 

0.1 

1 

315 

306 

—2 

0.3 

-  3 

0.8 

-1 

0.3 

360 

360 

-4-1 

0.5 

7 

—  3 

1.1 

31 

-1 

0.5 

16 

405 

414 

-4  2 

1.3 

—  2 

0.6 

0 

0.3 

450 

468 

-4-3 

1.6 

56 

0 

0.6 

11 

0 

0.5 

4 

495 

52? 

-4-1 

0.9 

2 

1.5 

1 

0.6 

540 

516 

—2 

0.7 

16 

2 

1.8 

77 

0 

0.8 

53 

585 

6S0 

-3 

2.1 

2 

1.0 

0 

0.5 

630 

684 

-6 

2.6 

120 

-  1 

0.8 

24 

0 

0.3 

8 

675 

738 

-4 

1.4 

-  6 

2.5 

-2 

1.1 

720 

792 

0 

1.1 

32 

—  8 

3.0 

147 

—4 

1.3 

63 

765 

846 

5 

3.4 

-  4 

1.6 

—2 

0.8 

810 

SOO 

5 

4.0 

207 

0 

1.2 

42 

0 

0.4 

.  13 

855 

954 

2 

2.2 

4 

3.8 

2 

1.5 

900 

1008 

-1 

1.4 

51 

4 

4.3 

239 

2 

1.9 

103 

945 

1066 

-5 

4.6 

2 

0.7 

1 

0.2 

9S0 

1120 

—8 

5.4 

319 

~  2 

1.5 

65 

-1 

0.5 

20 

1035 

1174 

-5 

3.0 

-10 

5.0 

—4 

2.0 

1080 

1228 

-4-4 

4.5 

76 

-10 

5.0 

353 

.  _4 

2.5 

152 

1125 

1212 

-4  8 

6.1 

—  6 

3.6 

—3 

2.2 

Saturn 

225 

120 

1  1 

0.2 

0 

0.4 

0 

0.2 

270 

256 

4-  3 

0.6 

0.6 

-  1 

1.0 

16 

0 

0.4 

6 

315 

392 

-4-  9 

1.4 

-  1 

1.1 

-1 

0.5 

360 

528 

-4-19 

1.4 

72 

4  6 

0.9 

4 

-4-1 

0.4 

1.4 

405 

664 

^16 

1.2 

-1 20 

2.4 

4-7 

1.0 

450 

800 

4  1 

1.2 

6.4 

-4  22 

3.4 

154 

4-9 

1.4 

64 

495 

936 

-  9 

4.0 

4 12 

2.6 

45 

1.2 

540 

1072 

-17 

6.2 

299 

—  9 

1.4 

15 

-4 

0.6. 

6 

Uranus 


225 

213 

4-  2 

3.3 

7 

-  1 

4.1 

8 

-  1 

2.0 

270 

587 

4-14 

2.5 

60 

—  7 

1.5 

12 

-  3 

1.0 

315 

910 

4  31 

52 

188 

-23 

4.8 

150 

-11 

4.8 

where  by  p  we  denote  the  error  in  a  single  step  due  to  the  neglf'cted 
disturbances,  and  otO  has  the  same  meaning  as  in  (l)  . 

In  the  computation  by  Formulas  (4)  for  all  three  planets  pis 
taken  as  the  maximum  disturbance  from  Neptune  and  Pluto;  that  is, 

I  •  10“*,  which,  of  course,  gives  a  strong  overestimate  for  Jupiter 
and  Hatum  (Table  2,  columns  5,  8,  ll) . 
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?>.  'A'e  shall  give  one  more  example.  Let  us  try  to  estimate  the 
error  due  to  round-off  in  the  eoordinates  of  the  larger  planets 
published  in  "Astronomical  Papers"  in  1951.  We  shall  make  the  esti¬ 
mate  by  the  crude  formula  obtained  from  (l)  under  the  follovdng 
assumptions:  we  neglect  and  E^.  in  comparison  with  we 

assume  equal  to  1.  Then  (l)  will  have  the  form 

.  (5) 

If  (j  is  taken  to  be  I  •  10~‘\  then  after  1000  steps  of  integration, 
which  corresponds  to  an  interval  of  time  greater  than  100  years  in 
this  case,  the  error  in  the  coordinates  of  the  planets  will  be 
approximately  equal  to  1 •  10~®  ;  that  is,  the  published  coordinates 
of  the  larger  planets  are  free  from  round-off  error. 

Thus,  the  examples  cited  above  show  that  the  formulas  derived  by 
B.  P.  Myanhin  for  the  estimate  of  errors  due  to  round-off  errors  is 
entirely  suitable  for  practical  use. 

Tile  estimate  (l)  reflects  the  oscillatory  character  of  the  error 
and  gives  a  comparatively  .smaller  overestimate  (as  a  rule,  by  less 
than  a  factor  of  10).  Furthermore,  it  shows  that  after  1000  steps  of 
integration  no  more  than  five  digits  are  lost  in  the  sought  for 
values  due  to  round-off  errors. 

But  as  for  the  error  created  by  the  neglected  disturbances 
(Tabic  '2),  the  estimate  (4)  whi  ..h  was  used  for  them  must  be  con¬ 
sidered  unsatisfactory,  since  it  docs  not  take  into  account  the 
periodic  character  of  the  disturbances.  This  estimate  gives  a  practi¬ 
cally  acceptable  result  only  if  during  the  entire  integration  on  a 
large  part  of  it  the  disturbances  are  constant. 
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